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The analytic equation of state of nonideal Coulomb plas- 
mas consisting of pointlike ions immersed in a polarizable 
electron background [G. Chabrier and A. Y. Potekhin, Phys. 
Rev. E 58, 4941 (1998)] is improved, and its applicability 
range is considerably extended. First, the fit of the electron 
screening contribution in the free energy of the Coulomb liq- 
uid is refined at high densities where the electrons are rela- 
tivistic. Second, we calculate the screening contribution for 
the Coulomb solid (bec and fee) and derive an analytic fitting 
expression. Third, we propose a simple approximation to the 
internal and free energy of the liquid one-component plasma 
of ions, accurate within the numerical errors of the most re- 
cent Monte Carlo simulations. We obtain an updated value 
of the coupling parameter at the solid-liquid phase transition 
for the one-component plasma: r m = 175.0 ± 0.4 (1<t). 

PACS numbers: 52.25.Ub, 52.25.Kn, 05.70.Ce, 97.20.Rp 



I. INTRODUCTION 

Fully ionized electron- ion plasmas (EIP) are encoun- 
tered in laboratory experiments, in stellar and planetary 
interiors, in supernova explosions, etc. From the theoret- 
ical point of view, the free energy of fully ionized EIP pro- 
vides the reference system for models aimed at describing 
the thermodynamic properties of partially ionized plas- 
mas. Thus the studies of EIP are of both theoretical and 
practical interest. 

In a previous paper p] we have calculated thermody- 
namic quantities of Coulomb plasmas consisting of point- 
like ions immersed in a compressible, polarizable elec- 
tron background and devised analytic fitting formulae 
for these quantities. The calculations were based on a 
linear-response theory for the ion-electron (ie) interac- 
tion, which is valid as long as the typical ie interaction 
energy (Ze) 2 /2ao (where ao is the Bohr radius, and Ze 
is the ion charge) is smaller than the kinetic energy of 
the electrons. This condition is fulfilled either at tem- 
peratures T > 10 5 Z 2 K or at densities p > AZ 2 gem -3 , 
where A is the ion mass number. For the nonrelativis- 
tic regime, i.e., at densities p <§; 10 6 gem -3 , finite- 
temperature effects were included in the electronic di- 
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electric function, as well as the local-field correction aris- 
ing from electron correlation effects, following the model 
developed in Ref. H. In the relativistic regime, similar 
calculations were done using the Jancovici H dielectric 
function. 

Since the electron screening is weak at high densities, 
and since the bulk of calculations have been performed 
using the nonrelativistic model, our fit for the ie contribu- 
tion was not very accurate at p > 10 6 gem -3 , where the 
electrons are relativistic. Because of the same weakness 
of the screening, this inaccuracy in the ie contribution 
at high p did not deteriorate the overall accuracy for the 
excess part of the free energy, which sums up the ion-ion 
(m), electron-electron (ee), and ie contributions. There 
are, however, physical problems which require an accu- 
rate evaluation of the ie part even at high densities (an 
example is mentioned below). In this paper we present 
a modification of the analytic formula ffl for the ie free 
energy which improves significantly the accuracy in the 
domain of relativistic electrons, keeping unchanged the 
previous nonrelativistic results. 

Second, we calculate the ie part of the free energy 
for a Coulomb solid, where the ions form either a body- 
centered-cubic (bee) or face-centered-cubic (fee) lattice. 
The calculation is performed in a perturbation approxi- 
mation, which is accurate because the screening is weak. 
We employ an analytic expression for the ion structure 
factor S(k) of a Coulomb crystal, obtained in Ref. || in 
the harmonic approximation for large wave numbers k 
outside the first Brillouin zone. For small k, we supple- 
ment it by an exact limiting form of S(k). We evalu- 
ate the screening contribution for both the classical and 
quantum harmonic crystals and construct a fitting for- 
mula which accurately reproduces our numerical results. 

The above mentioned improvements of the equation 
of state are significant at densities p > 10 6 gem -3 . 
Such densities cannot be reached in the laboratory, but 
they are commonly encountered in the interiors of white 
dwarfs and envelopes of neutron stars (e.g., Ref. fa]). 

In addition, we present simple formulae for the excess 
internal and free energies of a classical one-component 
plasma (OCP) liquid, which take into account the most 
recent Monte Carlo (MC) results H0, and which are ac- 
curate for any values of the Coulomb coupling parameter 
from the gaseous phase to the dense liquid regime. An- 
alyzing various results for the free energy of the OCP 
liquid and solid, we revise the value of the coupling pa- 
rameter at the solid-liquid phase transition. 



In the next section, we describe the basic parameters 
of the EIP. In Sec. |l|, we consider the OCP liquid and 
determine its freezing point. In Sec. IV, we present an 



improved fit to the free-energy contribution due to the 
electron screening in a Coulomb liquid. In Sec. M, we 
evaluate an analogous contribution for a Coulomb solid 
and fit it by an analytic expression. 



II. THERMODYNAMIC PARAMETERS 

We consider EIP consisting of pointlike ions and elec- 
trons. The basic dimensionless parameters are the elec- 
tron density parameter r s and the ion coupling parameter 
Y: 
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Here, ks is the Boltzmann constant, a e = (^iTUe)" 1 ' 3 
is the mean inter-electron distance, a — (^irrii)~ 1 ' 3 = 

a e Z 1 ' 3 is the mean inter-ion distance, and n e (rii) de- 
notes the electron (ion) number density. T e has a mean- 
ing of the coupling parameter for nondegenerate elec- 
trons. 

Quantization of the ionic motion is important if T <C 
T p = fuVp/ks, where ui v — (4irZ 2 e 2 rii/mi) 1 > 2 is the ion 
plasma frequency, rrii being the ion mass. A correspond- 
ing dimensionless parameter is 
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= T p /T = r y/3/Rs, 



where 
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is the ion density parameter. We neglect ion quantum- 
exchange effects, which is justified if Rs ^> T (see, e.g., 
Ref. |). 

The electrons are characterized by the degeneracy pa- 
rameter 9 and the relativity parameter x r , 



= T/T F 



x r =PF/(m e c), 



(4) 



where Tp is the Fermi temperature, c is the speed of 
light, and pp = Ti (37r 2 n e ) 1 ' 3 is the Fermi momentum. 
The electron screening properties are determined by the 
Thomas-Fermi wave number 
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where \x is the electron chemical potential. 

For these parameters, the following estimates are ac- 
curate within 0.005%: 

x r w 0.014005 r^ 1 « 1.0088 (p 6 Z/A) 1/3 , (6) 

22.547 x r „ i 5930 
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where pa = p/(10 6 gem' 3 ) and Ta = T/(10 6 K). In the 
nonrelativistic plasma (x r <C 1), 6 ~ 0.543 r s /T e . In 
the ultrarelativistic case (x r ~^> 1), 9 ~ (263 r e ) -1 . If 
the electrons are nondegenerate (9>1), fexFOe ~ V3r e . 
For strongly degenerate electrons (8 -C 1), 



fc TF a e « 0.185 (l + x~ 2 ) 1/4 . 



(8) 



The ion quantum parameter r\ is expressed through x r 
and r as 



7] w 0.3428 r y/x~r~Z- 7/6 A- 1/2 . 
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Within the aforementioned approximation of weak 
electron-ion coupling, the total Hclmholtz free energy 
.Ftot can be written as 
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where F^ d ' denote the ideal free energy of ions and elec- 
trons, respectively, and the last three terms represent 
an excess free energy arising from interactions. F-£ is 
the free energy of an ideal Boltzmann gas. For the elec- 

trons at arbitrary degeneracy and relativism, F id can 
be expressed through Fermi-Dirac integrals and approx- 
imated by analytic formulae fy]. An analytic parameter- 
ization for the nonideal (exchange and correlation) part 
of the free energy of the nonrelativistic electrons, F" e r , 
has been given in Ref. ||. For the relativistic electrons, 
the exchange free energy F*° l has been given, e.g., in Ref. 
fllPf , while the correlation corrections are negligible be- 
cause they contain an additional small factor ~ af In |a/| 
||ll|| , where a/ rs 1/137 is the fine-structure constant. In 
practice, we use the following interpolation between the 
nonrelativistic and relativistic regimes: if r e > 0.07 and 
r s < 0.13, we set 



F,„ 
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{l-Z)FZ + ZFf, 
exp[-(r e /0.07-0.9)" 2 
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(0.13/r s - 0.9)" 



otherwise we set F ee = F™. The interpolation is suffi- 
ciently smooth, because F" e r and F^ el closely match each 
other at the chosen boundary between the two regimes. 
In the following sections we consider the last two terms 
in Eq. (fl0|), which represent the excess free energy of an 
OCP of ions and the contribution due to the ion-electron 
interactions, respectively. 



III. OCP AND MELTING TRANSITION 

Liquid and solid phases of the OCP have been studied 
extensively by various analytic and numerical methods. 
All the thermodynamic functions of the classical OCP 
can be expressed as functions of the only parameter T. 
At f < 1, a diagrammatic cluster expansion yields 



u,, 



V3 



j^3/2 QF 1 ^ 



Nik B T 2 

-r 9/2 (1.6875 V3 lnT - 0.23511 



ln(3r) 



Cg _ 1 
2 3 
(12) 



where Cg = 0.57721 ... is the Euler constant. Here, the 
first term is the Debye-Hiickel energy, the second one is 
due to Abe p2| , and the ~ T 9 / 2 one is due to Cohen 
and Murphy JL3| . Since Fa vanishes at high T, it can be 
obtained from Uu by integration: 
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The above analytic expansion is not applicable for 
r > 1. The most accurate up to date numerical results 
for the internal energy of the liquid OCP at 1 < V < 200 
have been obtained by MC simulations by DeWitt and 
Slattery M and by Caillol M. These authors have also 
constructed analytic fits to their data with the standard 
deviations comparable to the numerical MC noise. Un- 
fortunately, these fits cannot be extended to small T, 
which hampers obtaining the free energy by Eq. (13). 
On the other hand, the hypernetted-chain (HNC) result 
for Fa at r = 1 is slightly inaccurate because the HNC 
approximation neglects the so called bridge functions in 
the diagrammatic representation of the interactions. To 
circumvent the difficulty, DeWitt and Slattery |Eij used 
small differences between HNC and MC at T = 0.8 and 
0.6 to get the corrected value of fu(T = 1) = —0.4368. 

We propose a different approach. We consider the pa- 
rameterization 
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where A3 = — v3/2 — A\j\fA^,. The terms in square 
brackets have been used in Ref. El, the term with B\ 
provides an adjustment of the fit to the MC data at 
large T, and the last term adjusts to Eq. ( |12|) at small V. 
The best-fit parameters with respect to the data BR] are 
given in Table |l|. Then the free energy can be obtained 
fromEq. (ph: 
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The corresponding expression for heat capacity is 
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FIG. 1. Upper panel: comparison of the fit (113) 

(solid line) with the Debye-Hiickel (DH), Abe Q, and Co- 
hen-Murphy |13| (CM) approximations (dot-dashed lines), 
with the MC results (circles) and the fit (dotted line) of 
Ref. [pi (DWS), and with some of our HNC results (trian- 
gles). Lower panel: residual differences between the fit (113) 
and (i) the analytic expansion ( |l2| ) (dot-dashed line), (ii) 
results of HNC calculations (triangles), (iii) MC results of 
Ref. Id] (open circles), and (iv) numerical results of Ref. M 
(MC+extrapolation) . 



Comparison of Eq. (15) with Eq. ( |12] ) at T < 1 and 
with the MC data from Refs. MM at T > 1, supplemented 
by some of our HNC calculations, is given in Fig. [j]. The 
upper panel displays the ratio uu/T 3 ' 2 (which is constant 
in the Debye-Hiickel approximation). The magnitude of 
the possible error is demonstrated by the lower panel. 
Here, the dot-dashed line shows the difference between 
the approximation ( |15|) w ith the second set of param- 
eters and expansion (|12[), while various symbols show 
residual differences between the same approximation and 
numerical (HNC and MC) results. The distribution of 
the residuals around zero looks irregular, which indicates 



that they represent a numerical noise of the MC calcula- 
tions rather than an error of the fit (y_5|) . In addition, we 
have checked that the difference between our fit to the 
free energy, Eq. (|16|), and the one in Ref. H (at T > 1) is 
of the order of the aforementioned small uncertainty in 
/«(T = 1). 

More complicated interpolations between the low- and 
high-r limits were proposed previousl y ||15| , |lq| . By con- 
struction, they reproduce exactly Eq. (|12|) at Y — > and 
the fits to MC results at T ^ 1. Compared with the 
present fit, however, those interpolations have somewhat 
larger differences from the HNC results at 0.1 < V < 1. 
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FIG. 2. Difference between the free energy of the solid 
OCP given by a 3-parameter fit of Ref. |2G] and parameter- 
izations for the liquid OCP according to Refs. |19] (SDD; 
short-dashed line), B (dotted line), Eq. ( Jlq ) (solid line), and 
for the solid according to Dubin |l£J] (dot-dashed line). The 
long-dashed lines marked "FH" correspond to the 4-parame- 
ter fit and to the ±lc- uncertainty of the 3-parameter fit in 
Ref. @. 

The freezing of Coulomb OCP liquid into a bcc crystal 
occurs when the free energy of the solid becomes lower 
than that of the liquid at T = r m . Nagara et al. Jl7[ ] 
and Dubin |lq | , having improved a previous treatment of 
anharmonic corrections to the free energy of the Coulomb 
crystal, obtained Y m = 172 ± 1. However, these authors 
employed an older fit [[19| (SDD) for the liquid. Figure g 
shows the differences between fa for the solid and liquid 
OCP given by various parameterizations. For the solid, 
we have adopted the three-parameter fit by Farouki and 
Hamaguchi |£0[ to their molecular dynamics simulations 
in the range 170 < Y < 400. The horizontal long-dashed 
lines correspond to the standard deviation of that fit. 
The line between them represents a four-parameter fit 
P0j| in the same Y interval. The dot-dashed line shows 
the difference between the fit of Ref. pi and that by 
Dubin Q . The value of Y m indicated above is given by 
the intersection of the latter line with the short-dashed 
one (SDD). Using updated results for the OCP liquid 
(either Ref. || or our Eq. (|16|) , represented by the dotted 
and solid line, respectively) and the OCP solid pQ] , we 
obtain Y m = 175.0 ±0.4. 



IV. ELECTRON SCREENING IN A COULOMB 
LIQUID 

We now consider electron polarization effects in the 
EIP. In the previous paper Q, we have calculated i*i e 
using the model developed in Ref. |g] for nonrelativistic 
EIP. The HNC equations have been solved numerically 
for an effective screened inter- ion potential V d q, which is 
the sum of the bare ionic potential and the induced po- 
larization potential, to obtain Fu+Fi e and corresponding 
contributions to the internal energy (Ua + Ui e ) and pres- 
sure (Pa + Pie). The same equations solved for the bare 
Coulomb potential give Fa, Uu, and Pa. The difference 
represents the screening (ie) part. Inclusion of the finite- 
temperature effects in V c s provides a correct treatment 
of the thermodynamic quantities over a wide range of 
values of Y from the Debye-Hiickel limit T <C 1 to the 
strong-coupling limit T ^> 1 for various r s and Z . 

Relativistic calculations have been performed employ- 
ing the same HNC technique but with the Jancovici 0] 
dielectric function e(k, x r ), which is appropriate at strong 
degeneracy (8 <C 1) and arbitrary x r . The results are in 
good agreement with those obtained by Yakovlev and 
Shalybkov pi[ , who have used an equation 

/* = — f- = / S(k,Y)[e(k,x r )-l]dk, (18) 

Nik B T 7T J 

where S(k,Y) is the static structure factor of ions (i.e., 
the Fourier transform of the ion radial distribution 
function). Equation (18) has been derived by Galam 
and Hansen pl[ using a thermodynamic perturbation 
scheme, which can be represented as an expansion in 
powers of /ctf- We have repeated the calculations |ll| ] 
using a more recent and accurate S(k, Y) ||22f than in the 
original work; the change in /j e due to this update does 
not exceed 4%. 

Note that Eq. (n8h differs from the standard first-order 
perturbation approximation by a replacement of [1 — 1/e] 
by [e — 1]. The resulting difference « (kTpa) 3 /6 has the 
same order of magnitude as the second-order perturba- 
tion correction pl| . Our HNC calculations with the Jan- 
covici dielectric function at x r % 1 and T > 1 coincide 
within 2% with the results of Ref. JO], whereas the sub- 
stitution of [1 — 1/e] in Eq. ( |18| ) yields a considerable 
difference: for example, for Y = 1 and Z = 26 this differ- 
ence approaches 40% even at large x r . We conclude that 
the approximation (|18| ) is very accurate at high densities. 

The screening contribution to the free energy of the 
Coulomb liquid at < r s < 1, < T < 300, and 1 < 
Z < 26 has been fitted by the expression |p| 
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where c DH = (Z/y/3) [(1 + Zf' 2 - 1 

act transition to the Debye-Hiickel limit at T — ► 0, 

ctf = (18/175) (12/^) 2 / 3 Z 7 / 3 (1 - Z-V3 + 0.2 Z-V2) 



reproduces the Thomas-Fermi limit j23| at Z — > oo, the 
parameters a = 1.11 Z 0A75 , b = 0.2 + 0.078 (InZ) 2 , and 
z/ = 1.16 + 0.08 In Z provide a low-order approximation 
to Fi, e for intermediate r. and T, and the functions 
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improve the fit at relatively large r s . The results of our 
nonrclativistic finite-temperature HNC calculations are 
reproduced by setting hi and hi equal to unity; these 
factors come into play in the relativistic case. 
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FIG. 3. Calculated (filled circles) and fitted (solid lines) 
normalized contribution to the internal energy due to polar- 
ization, Ui e , as function of V at different values of x r , Z = 6 
(left panel) and as function of x r at two values of Z, T — 150 
(right panel]. For comparison, approximations M (dotted 
lines) and |llj (dashed lines) are also shown. 

In the latter case, the asymptotic behavior of Eq. ( |l9| ) 
at r — ► oo should change from /j e oc Tr s to /j e oc 
Tr s y/1 + x 2 . This is achieved simply by setting hi = 

(1 + x| 

limit |2 
actly. 

The factor hi is devised to correct the fit at finite Z 
in the relativistic domain. A form chosen previously M 
was not very accurate, as illustrated by the dotted lines 
in Fig. for the internal energy 



1 ' 2 . Then the zero-temperature Thomas-Fermi 
j] (r» < 1, f -> oo, Z — > oo) is reproduced ex- 



A more accurate relativistic correction reads 

1 + 3:2/5 



hi(x r ) 



1 + 0.18 Z-y^Xr + 0.37 Z-V*x* + xljh ' 



(21) 



The resulting Ui e [Eq. (pp|)] is plotted by the solid lines in 
Fig. |3j. There is now a good agreement with the thermo- 
dynamic perturbation expansion pd| ] at large T for any 
x r , without deteriorating the accuracy of the old fit in the 
nonrclativistic domain. Quantitatively, for 1 < T < 100 
and x r < 0.25, the difference between the fit and the 
HNC results is typically 2-3%, with maximum 8% for 
Z = 1, r = 100 and r s = 2.074 (the maximum r s value 
used in the calculations). Note that the model of EIP has 
only marginal physical relevance at such large values of r s 
and r because of the incipient bound-state formation. On 
the other hand, at very strong coupling (r > 100) and 
relativistic densities (x r > 0.1), the results of Ref. pl[] 
and of our relativistic HNC calculations are reproduced 
by our fit with typical deviation of 1-3% (maximum 4.3% 
at Z = 6, T = 100, and x r = 10). 
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FIG. 4. Absolute values of the heat capacity of fully 
ionized liquid carbon at T — 10 K and 10 K. Dotted 
curves show the contributions of the electrons (heavy line - 
ideal Fermi-gas contribution, light line - exchange and cor- 
relation correction), dashed lines - contributions of the ions 
(long dashes - ideal-gas part, short dashes - correlation part), 
and dot-dashed curves - ion-electron (polarization) correc- 
tion. The latter curves end at V = 175. The dips on the ee 
and ie curves signify a change of sign. For the ii and ie contri- 
butions, present approximations (heavy lines) are compared 
with those in Refs. § (DWS) and || (YS) (light lines).The 
heavy solid line shows the sum of all terms. 

The heat capacity per ion in units of ks, Cy /Niks, 
of the classical EIP liquid is shown in Fig. [| for Z = 6, 
T = 10 6 K and 10 7 K. These plasma conditions can oc- 
cur, for example, in interiors of some giant stars or in 
accreted envelopes of neutron stars J24j. Various con- 
tributions, shown in the figure, correspond to separate 
terms in Eq. (|l0|). At relatively low densities, the main 
contribution is that of the electrons, with the limiting 



value \Z = 9. With increasing p, the electron gas be- 
comes degenerate, and its heat capacity decreases. Then 
Cy is determined by the ion liquid. The Coulomb (it) 
contribution slightly exceeds the kinetic one (I) near 
freezing. According to the equipartition theorem, in a 
classical ionic crystal the potential and kinetic contribu- 
tions are each equal to | (apart from small anharmonic 
corrections) . This means that freezing is accompanied by 
a drop of Cy, equal to the excess of the potential contri- 
bution over the kinetic one in the ionic liquid just before 
freezing. We see however that this excess (and hence the 
drop) is not large. 

The values of Cv,u determined by Eq. ( p"7| ) (thick 
dashes) and derived from the fit in Ref. M (thin dashes) 
are close to each other near the freezing. With decreas- 
ing density, however, a large difference develops, which 
is natural because the formula in Ref. [|6| is not applica- 
ble at small T. Of the same origin is the striking dis- 
crepancy between the approximations for Cy,ie derived 
from Eq . dl3 ) (thick dot-dashed curve) and from the fit 
in Ref. |11| (YS), seen at low p. In this domain, our 
fit describes the change of sign of Cy.ie from negative 
in the strong-coupling regime to positive in the Debye- 
Hiickel domain. However, an appreciable difference with 
Ref. O] persists even at large p, where both fits describe 
fi e equally well (within uncertainties in the structure fac- 
tor). This reflects insufficient accuracy of the present-day 
determination of the functional form of S(k, T) for the 
strongly coupled Coulomb liquid. 



V. ELECTRON SCREENING IN A COULOMB 
SOLID 



ment would consist in calculating the dynamical matrix 
and solving a corresponding dispersion relation for the 
phonon spectrum. The first-order perturbation approxi- 
mation for the dynamical matrix of a classical Coulomb 
solid with the polarization corrections, based on an ef- 
fective inter-ion potential, was derived by Pollock and 
Hansen |26|. In a quantum crystal, strictly speaking, one 
would have to consider the electron-phonon interactions, 
in order to calculate the perturbed spectrum. 

As mentioned above, the polarization of the electron 
gas is weak at the high densities we are interested in. This 
suggests a simpler, semiclassical perturbation approach 
to evaluate the polarization corrections. The ionic crys- 
tal without ie interactions is a natural reference model. 
Note that the effective inter-ion potential in the adia- 
batic perturbation approximation 26 is just the electro- 
static potential, common to the liquid and solid phases. 
The difference of this potential from the bare Coulomb 
potential can be considered as perturbation. Then we 
can apply the Galam-Hansen |2l| perturbation theory, 
which is based on the exact expression for the free en- 
ergy involving an integration over a coupling parameter 
related to the "strength" of the perturbation. Thus we 
recover Eq. ( |J8| ) in the case of solid, with S(k) replaced 
by (47r) _1 J Sjk)d£l, where dQ is a solid angle element in 
the direction of k. 

The resulting polarization correction ( fig ) does not take 
into account quantum aspects of the ie (electron-phonon) 
interactions, but it allows us to study effects arising 
from quantum modifications of the ion-ion correlations. 
These correlations are described by the structure factor 
S, which depends in this case on k, T, and -q. 



A. Perturbation approximation 

At high densities and below a certain temperature, the 
ionic Coulomb plasma forms a Wigner crystal. For ex- 
ample, interiors of cool white dwarfs fl25fl are expected 
to be in the solid state. The cooling is governed essen- 
tially by the compressibility and heat capacity of their 
interiors, whose central regions are compressed to rel- 
ativistic densities. In that case, the main contributions 
to the internal energy (the zero-temperature electron-gas 
kinetic energy and the ion electrostatic part) do not de- 
pend on temperature, so that the heat capacity is entirely 
determined by small temperature-dependent corrections. 
Therefore, evaluation of the polarization corrections for 
the Coulomb solid is important for astrophysical appli- 
cations. 

Since the maximum ion frequency in the solid is quite 
small compared to the electron plasma frequency, one 
can use the adiabatic (i.e., Born-Oppenheimer) approx- 
imation, which allows to decouple the electron and ion 
dynamics. Even so, a calculation of the thermodynamic 
functions of a Coulomb solid with allowance for the ie 
interactions is a complex problem. A rigorous treat- 



B. Structure factor 



In a crystal, the static structure factor is given by 
S(k,T,n) = ij e ik '( R, - H *)(e il "V ik, '') r , (22) 

where vn is an operator of ion displacement from an equi- 
librium lattice position R^, and (. . .)t denotes the canon- 
ical average. The structure factor (p2|) can be decom- 
posed into elastic (or static-lattice) and inelastic parts, 

SQi, r, v ) = s'(k, r, v ) + s"(k, r, v ). (23) 

The elastic part is (e.g., Ref. Q) 

S' sol (K T, rj) = (2nf ni e' 2W ^ r ^5(k G), (24) 

G 

where ^ G denotes a summation over all reciprocal lat- 
tice vectors G but G = 0, and e~ 2W = (exp(ik • u)) T is 
the Debye- Waller factor. In isotropic (e.g., cubic) crys- 
tals, one has 



2W(k,T,r 1 ) = r 2 T CT,r ] )k 2 /3 J 



(25) 



where r 2 -, 



(u 2 )t is the mean-squared ion displacement 



(cf. ]27|]). In a harmonic crystal, 



2 a27 ? 



M-i 



1 



uj v exp(huj L ,/kBT) — 1 



ph. 



. (26) 



where v = (q, s), s = 1,2,3 enumerates phonon modes, 
q is a phonon wave vector, u> v is the frequency, (. . .) p h 
denotes averaging over phonon wave vectors and polar- 
izations, and fx n = ((uj l ,/ujp) n )pi l . In the classical limit 
(rj — > 0), r\ = /i_2fl 2 /r; and in the quantum limit 
(77 — ► oo), r|, ~ /i_ia 2 ry/(2r). Numerical values of /n_i 
and /z_2 are given in Table O. At arbitrary r), a con- 
venient analytic approximation to r\ is provided by a 
model of the harmonic Coulomb crystal p8| which treats 
two acoustic modes as degenerate Debye modes with 
w u = aoj p q/qBz, where gez = (Gir^rii) 1 ' 3 is the equiv- 
alent radius of the Brillouin zone, and the longitudinal 
mode as an Einstein mode with w v = ^ujp. Accuracy of 
this model for the thermodynamics of the bcc Coulomb 
crystal has been demonstrated in Ref. [M , where the val- 
ues a = 0.399 and 7 = 0.899 have been derived from the 
requirement that the model should reproduce the exact 
values of /i_2 and /12 = |. For the fee lattice we obtain 
a = 0.413 and 7 = 0.892. Using this model, we can 
calculate the second term in Eq. (Eq), which yields 



a 

T 



V-iV 



V 



37 (ei"> - 1) a 3 r) J e 



ari tdt 



1 



(27) 



This approximation ensures the correct classical and 
quantum limits. Between these limits, the maximum de- 
viation from accurate numerical results |33| ] reaches 1.6% 
at 77 « 9 for both bcc and fee lattices. 

According to Eq. (23), Eq. (18) can be rewritten as 



lie — lie T /, 

f - 
J ie 



11 



(28) 



T-T! e{G tn:l - *M-2W(G, r, ,)], (29) 



G 



(GaY 



•p /*00 

/;: = -—/ S"(k,T,r,)[e(k,x r )-l]dk. 
k Jo 



(30) 



The inelastic part of the structure factor of a harmonic 
crystal reads Hi 



S" = 



-2W 



R 



ikR 



^v(R)k 2 



1 



cos(q • R) 



3ft / (k-e,) 2 

2?7ti \ k 2 lo v t&znhifuJv/lkBT) / h 



(31) 
(32) 



(where e„ is a phonon polarization vector). A straight- 
forward use of this expression is impractical because of a 
slow convergence of the sum. For this reason, we employ 
the approximation 01: 



S"(k,r,r)) w S '(fc,r,77) = l-exp[-2W(fe,r,77)]. (33) 

As argued in Ref. H], this approximation is good for use 
in integrals over k at k > qbz- In papers addressed to 
the transport properties of Coulomb plasmas @,Q, in- 
tegrals over k were truncated from below at k = gez- In 
Eq. (p0|), however, it is essential to recover the correct 
limiting behavior of S"(k) at k — * 0, since [e — 1] oc k~ 2 
becomes large in this limit. Therefore we use a piecewise 
approximation: 



S" = 



(k>h), 

si' = ifc 2 \d 2 s"/dk 2 ] k ^ (k < k x ), 



q/l 



(34) 



where the parameter k\ will be determined below. The 
exact result for classical Coulomb plasmas pl[ reads 
S'i(k) = (fca) 2 /(3r). In general case, S'i(k) can be found 
from Eq. (|3l|). At small k, the expression in the square 
brackets in Eq. ( J3l| ) can be replaced by v(R)k 2 , which 
corresponds to the one-phonon approximation. Chang- 
ing the order of averaging and summation, we see that 
the summation yields delta function i5(k±q— G); there- 
fore q = ±k as long as k < minG ~ 2^bz- Hence, only 
the longitudinal phonon mode contributes in this limit. 
The frequency of this mode in a Coulomb crystal tends 
to ujp at small q (e.g., ||26| ), which enables us to perform 
averaging in Eq. (B3). Finally we obtain 



£i(*,r,T7) 



(fcq) 2 

6r tanh(r//2) ' 



'/ 



(35) 



In order to test our approximation (|34|) and to find 
the optimum value of k\ , let us consider the electrostatic 
energy U c \. st of a Coulomb crystal, 



U, 



Wel-st 



cl-st 



= — / \S(k) - 1] dk = v! 



N,k B T 
where, according to Eqs. (|23|) and (|24|), 



3r 



G 



(Ga) 2 



is the static-lattice part. Baiko et al. 
that 



(36) 



(37) 



have shown 



u - = -c 
r 



M 



2a 2 



3 a 

7r 2tt 



where the terms not explicitly written are exponentially 
small at large I\ For the inelastic contribution, our model 



yields u" — u'q + ■ 



where 



To 

Ti" Jo 
To 

n Jo 
(ha) 3 



[£o'(fc)-l]dfc 



[S'{(k)-S' '(k)]dk 



3 Ta 

7T 2tt ' 



Ta 



18tt tanh(ry/2) 



fci- 



3tt r r T ki 
err- 



2r T 



V3 



(38) 



(39) 



On the other hand, in the harmonic lattice approxi- 
mation, £/d_ s t = —NiCM{Ze) 2 /a + U v /2, where the first 
term represents the energy of a perfect ionic lattice in 
uniform electron background, Cm being the Madclung 
constant (Table ph, and, from the virial theorem, the 
second term is one half of the vibrational energy of a 
harmonic crystal, 



U v = 3N l k B T 



>1 



uj p e v^ v /u p _ i 



P h 



2 



(40) 



We determine fci so as to recover the classical limit 
Ueist — — CmT 4- 3/2 at r\ = and T — ► oo. This yields 

1/3 

« 0.94. (41) 



<?BZ 



M-2 



M-2 - 1 



Figure || shows u i_ st calculated from Eqs. (p7|)-(p9|) for 
the bcc crystal at finite r\ and T (dot-dashed lines), com- 
pared with a calculation in which S" is set equal to Sq at 
any k (dotted lines) and with results of numerical calcula- 
tions Q| . We see that our modification of the structure 
factor at k < 0.94 qbz provides a significant improve- 
ment over the model without such modification (denoted 
as HL1 in Ref. ||). 
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FIG. 5. Normalized electrostatic energy of a bcc Coulomb 
crystal calculated using the approximate structure factor 
given by Eq. (B4J) (solid lines) compared with analogous calcu- 
lations with a model structure factor, Eq. (|33| ) (dotted lines), 
and with accurate numerical calculations J33J (triangles) . Up- 
per curve of each type or symbol corresponds to T — 200 
and lower one to F = 500. Long-dashed line displays the 
Madelung limit. 



C. Results 



Using Eqs. (|28|)-(|30|) and fl33|)-(|35|), we have calculated 
the polarization correction /j e over a wide range of pa- 
rameters: 80 < T < 3 x 10 4 , 10~ 2 < x r < 10 2 , and 



1 < Z < 92. Not all combinations of the considered pa- 
rameters are physically relevant; for instance, at Z = 1 
and large x r the ion-exchange effects neglected in our 
study become important. The use of such extended set 
of parameters, however, delivers robustness to a fitting 
formula presented below. 
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FIG. 6. Normalized polarization correction fi e to the 
free energy of a Coulomb crystal. Left panel: log 10 (— fi e /Y) 
against log 10 Y at Z — 6 for several values of x r . Right panel: 
log 10 (— fie/Y) against log 10 x r at Y — 10 3 for several values 
of Z. Solid lines: a classical solid; dot-dashed lines: quantum 
effects included. On the left panel, dotted lines show the re- 
sults for the classical solid with simplified structure factors, 
and dashed lines show the results in a liquid. 

Some results of our calculations are shown in Fig. [(j. 
Solid lines correspond to the piecewise approximation 
(|34| ) of the structure factor in the classical case (r\ — ► 0). 
Dashed lines on the left panel reproduce calculations in 
the liquid with S(k,T) from Ref. P2| . The upper and 
lower dotted lines at every value of x r show, respectively, 
the results of calculations with the inelastic part of the 
structure factor replaced by Sottas in the HL1 model of 
Ref. J32|) and by (as in Refs. P,pO|). Compared to these 
simplified approximations, the present model provides a 
smaller discontinuity of /j e at the freezing point (near 
the ends of the dashed lines). On the other hand, the 
divergence of the dotted curves towards smaller T shows 
that the result is still model-dependent. This model de- 
pendence disappears at T > 3000, since the static-lattice 
contribution becomes relatively large. 

In reality, at large values of T and small values of x r 
shown in Fig. pL the quantization of ionic vibrations be- 
comes important. This quantization considerably modi- 
fies the structure factor. This effect is taken into account 
by letting 77 to be finite in Eqs. (|27]) and (J35|). Results 



of the calculations, where r\ was determined from Eq. (Q) 
assuming A = 2Z, are plotted in Fig. g by dot-dashed 
lines. The curves on the left panel become flat as r\ be- 
comes large, which corresponds to an approximate pro- 
portionality fi e ex r. As a consequence, the polarization 
contribution to the specific heat, Cy,ie, goes to zero at 
large r\ (but remains one of the leading contributions, as 
shown below). 

The numerical results can be fitted by the expression 



fie = -foo(xr) r [i + A( Xr ) (Q(v)/ry] 



(42) 



where 



foo(x) = a TF Z 2 / 3 h ^/l + b 2 /x 2 , 
b 3 + a 3 x 2 



A(x) = 



1 + b 4 x 2 



Q{ri) = \/l + (OT) 2 , 

and parameters s and 61-64 depend on Z: 



1 + 0.01 (In Z) 3/2 + 0.097 Z~ 2 



b\ = 1 — a,\ Z 



-0.267 



0.27 Z~ 



1 



2.25 1 + a 2 Z 5 + 0.222 Z 6 
ZV3 



1 + 0.222 Z 6 



63 = 04/(1 + In Z), 

64 = 0.395 In Z + 0.347 Z~ 3/2 . 

The parameter otf, related to ctf in Eq. jlSJ), is cho- 
sen so as to reproduce the Thomas-Fermi limit [£3J at 
Z -> 00: a T F = (54/175) (12/tt) 1 / 3 a/ = 0.00352. The 
numerical parameters a \-a 4 and q are slightly different 
for bec and fee crystals; they are given in Table III. 

For a classical crystal, an average error of the fit is 1% 
for all Z, x r , and T, and the maximum error is 3.1% at 
Z = 92, r > 10 4 , and x r sw 2. In the quantum case 
(r) 7^ 0), the fit is accurate for Z > 3 only. In the range 
3 < Z < 30, an average error is 1%, and a maximum 3% 
occurs at Z = 3, T rs 100, and av w 2. 



D. Discussion 

The results presented in Fig. g indicate that, although 
the polarization in a Coulomb crystal is very weak, it 
does not vanish even at arbitrarily large T and x r . As in 
the case of strongly coupled liquid, fi e is roughly propor- 
tional to (/ctfo) 2 , which tends to a finite limit at rela- 
tivistic densities. The order of magnitude of the screening 
correction 6i e — Fi e /F to t for a classical Coulomb plasma 
at arbitrarily high densities is given by the Thomas- Fermi 
result |23[] , which is reproduced by Eq. ( |42] ) at T — > cxd 
and Z — >• 00: (Ji e « 0.004 Z 2 / 3 . Quantitative difference 
of the perturbation result at finite Z from the Thomas- 
Fermi limit is quite noticeable, ~ Z~ - 3 . 
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FIG. 7. Upper panel: difference between polarization cor- 
rections to the free energy in the solid and liquid phases at 
r = 175, x r = 1, 3, 10. Lower panel: Coulomb coupling 
parameter at the melting point when polarization corrections 
are taken into account. 



As mentioned in Sec. V A, our treatment of the screen- 
ing contribution is approximate. Nevertheless, we can 
use these results in order to demonstrate the importance 
of the polarization corrections. On the upper panel of 
Fig. 0, the difference Af ie between f ie values in the solid 
and liquid Coulomb plasmas at the OCP melting point 
r = 175 is plotted against Z for three values of x r . The 
largest x r = 10 represents virtually the ultrarelativis- 
tic limit. When compared to Fig. @, this plot shows that 
A/j e is sufficiently large to affect T rn . This effect is shown 
on the lower panel of Fig. M, where we have plotted our 
estimate of T m at x r = 1 and x r S> 1. Since A/j e re- 
mains finite at any x r , the classical OCP value r m = 175 
is never exactly recovered even at arbitrarily large p. 

Another important effect of the polarization correc- 
tions in the solid phase is that on the specific heat Cy- 
By differentiation of Eq. (B3) , we obtain 



u ie = -fooT [1 + A (1 - s/Q 2 ) {Q/V) 8 ] 



%£- '-* ( * 



(qr,f 



Q 3 



(43) 
(44) 



In a classical crystal, Cv,ie is only a small negative cor- 
rection to the total Cy rs 3Nikg. When T decreases 
much below T pi the heat capacity of an ionic crystal E9] 
goes to zero as Cy.i ~ l-6NikBTr 4 /(ari) 3 ex T 3 , whereas 
the ie contribution becomes positive and decreases as 



C 



V,ie 



Niksfoo sA (i? s /3 9 2 ) (1 - s)/2 ( W )- 1 ex T, (45) 



at the same rate as the heat capacity of a strongly de- 
generate electron gas |11 , 



Cy e ~ ZN t (k B T/m e c 2 ) ir 2 y/l + x 2 /x 2 r 



(46) 



Equation (|45|), derived from the fit J4^), agrees with 
the limiting expression at 77 — ► 00 which follows from 
Eqs. @, (H, (H, and (| 



C 



V,ie 



N,k r 



a 3 r) 



V e(G,x r ) - 1 



Lt- 3 




/•OO 


-1 


/ [e(fc,x r )- 


- 1] k 2 dk 


Jfci 


. 



(47) 



Thus Cy,ie becomes larger than Cv,i at sufficiently low 
T, which probably signifies that the thermodynamic per- 
turbation theory is violated at this T. 

The discussed effect is of anharmonic nature. Indeed, 
the harmonic approximation for the Hamiltonian leads 
to the Debye law Cy oc T 3 , regardless of inclusion of the 
polarization correction in the force matrix. Therefore 
the dependence Cy.ie oc T in Eqs. (^) and ( |47| ) is due 
to the use of the full Coulomb potential (not only its 
harmonic part) in the ie interaction energy, which has 
led to Eq. (|l|). 

It is also noteworthy that the modification of the OCP 
structure factor by the quantum effects renders Cy,ie pos- 
itive. A plain extrapolation of the ie contribution from 
the liquid regime into the solid would be completely in- 
appropriate, as it would result in a negative total heat 
capacity. 

The behavior of different contributions to the heat ca- 
pacity in the solid phase as function of p and T is il- 
lustrated in Fig. M Here we consider 12 C at 10 5 K and 
10 6 K. In the latter case (the bottom panel) one can see 
also the discontinuities at the liquid-solid phase transi- 
tion at p « 10 5 gem -3 , discussed above. As in a liquid, 
we can safely neglect the exchange correction, which at 
x r ^S> 1 is as small as — (a//27r) Cy. e ~ — 10 -3 CV. e . At 
relatively low densities, Cy is determined mainly by the 
ionic contribution. As T p becomes greater than T with 
increasing density, the phonon contribution to Cy freezes 
out rapidly, and Cy becomes determined by the degen- 
erate electron gas, polarized by the electric field of ions. 

This may have important consequences for astrophys- 
ical applications. In particular, the heat capacity of old 
white dwarfs, whose temperature is so low that their in- 
teriors are formed of quantum Coulomb crystals p5[ , 
may be substantially influenced by the polarization ef- 
fects IP. 
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FIG. 8. Absolute values of heat capacity of carbon at 
high densities for two values of T. The contributions of free 
electrons, ionic OCP, and electron-ion interaction are shown 
by dotted, thin solid, and dot-dashed lines, respectively. The 
thick curve shows the total value. The dashed part of the 
latter curve corresponds to the region where the thermody- 
namic perturbation theory used in the calculation of Cv,ie is 
not reliable. 

An improvement of the ii part enables us to deter- 
mine accurately the classical OCP melting point. The 
Coulomb coupling parameter at the phase transition is 
found to be T m rs 175, slightly larger than a previ- 
ously determined value. An improvement of the ie part 
in the liquid phase yields a better precision at densities 
p > 10 6 gem -3 , where the electrons are relativistic. Fi- 
nally, our estimates of the ie part of the free energy of 
a Coulomb crystal show that it is important for applica- 
tions. For example, our results demonstrate that it af- 
fects the melting of a classical Coulomb crystal and may 
contribute appreciably to the heat capacity of a quantum 
crystal. Since our calculations for the Coulomb solid are 
based on an approximate method and performed using 
an approximate structure factor, the latter results can 
be considered as estimates only. These estimates show, 
however, that the polarization corrections in Coulomb 
crystals are not as unimportant as it was often believed; 
they deserve to be studied further using more elaborate 
methods. 



VI. SUMMARY 

We have improved analytic approximations jl| for the 
contributions to the free energy of a Coulomb liquid due 
to the ii and ie correlations. In addition, we have sug- 
gested an approximation for the ie contribution to the 
free energy of a Coulomb crystal. 
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TABLE I. Parameters of Eq. @. 



data 



-Ax 



M 



Si 



B, 



-B 3 



Bi 



Rcf. 
Ref. 



II 0.9070 0.62954 4.56[-3] 211.6 1.0[-4] 4.62[-3] 
0.907347 0.62849 4.50[-3] 170.0 8.4[-5] 3.70[-3] 



a Powers of 10 are given in square brackets. 



TABLE II. Parameters of Coulomb crystals Q. 


lattice type fi-2 /i-i /ii 


Cm 


bcc 12.973 2.798 55 0.5113875 
fee 12.143 2.719 82 0.5131940 


0.895 929 255 682 
0.895 873 615195 


TABLE III. Parameters of Eq. 


(§)• 


lattice type at 02 0,3 


a i q 


bcc 1.1866 0.684 17.9 
fee 1.1857 0.663 17.1 


41.5 0.205 
40.0 0.212 
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